In [6]:
import skgeom as sg
import matplotlib.pyplot as plt
import numpy as np
import logging
import heapq
from euclid3 import *
from itertools import *
from collections import namedtuple
In [7]:
# _window function takes a list, lst and returns a zipped format of iterables prevs, items and next
def _window(lst):
prevs, items, nexts = tee(lst, 3) # Return n independent iterators from a single iterable, which is the list lst
prevs = islice(cycle(prevs), len(lst) - 1, None) # cycle gives an infinite cycle, for example lst = ['A', 'B', 'C'] cycle(lst) will give 'A', 'B', 'C', 'A', 'B', 'C', ...
#for prevs, islice slices from last element till infinity
nexts = islice(cycle(nexts), 1, None) #islice starts from index 1 till infinity
return zip(prevs, items, nexts) # returns zipped list, zip uses the list with shortest length (item) to do the pairing so for every item there exists previous and next
In [8]:
# A demonstration of how the itertools islice and cycle works in _window function
#vert_lst1 = [(100, 50), (200, 200), (250,40), (635,10), (635,250), (500,235), (625,425), (40,520), (40,50)] # Vasan Convention
poly_vert = [(40, 50), (40, 520), (625,425), (500,325), (635,250), (635,10), (250,40), (200,200), (100,50)] # Polyskel Convention
#vert_lst = [[100, 50], [200, 200], [250,40], [635,10], [635,250], [500,235], [625,425], [40,520], [40,50]]
iter_it=_window(poly_vert)
# the iterables can be iterated through just one time
for prev, point, next in iter_it:
print(prev,point,next)
In [9]:
# Sometimes we might get polygon vertices which lie on the same line or the same point might be repeated again. To avoid that, we normalize the contour. We also convert the vertices to a point instance of Euclid package for geometric manipulation
def _normalize_contour(contour):
contour = [Point2(float(x), float(y)) for (x, y) in contour]
return [point for prev, point, next in _window(contour) if not (point == next or (point-prev).normalized() == (next - point).normalized())]
In [10]:
norm_vert = _normalize_contour(poly_vert)
In [11]:
norm_vert
Out[11]:
In [12]:
#defining the cross product (AxB) for vectors A(ax,ay) and B(bx,by)
def _cross(a, b):
res = a.x * b.y - b.x * a.y
return res
In [13]:
# To avoid rounding off errors, we have two functions _approximately_equals(for vectors) and _approximately_same(for points). The function will return both the values are equal if the absolute difference between them are less than 1/1000th of the maximum value
def _approximately_equals(a, b):
return a == b or (abs(a - b) <= max(abs(a), abs(b)) * 0.001)
def _approximately_same(point_a, point_b):
return _approximately_equals(point_a.x, point_b.x) and _approximately_equals(point_a.y, point_b.y)
In [14]:
# __slots__ = () is used to avoid the creation of dictionary and faster access to attributes, it is mainly used along with namedtuple ( a built in datatype from collections packages), namedtuple allows us to access the tuple by both keyword and index
# __lt__(self,other), is a method of writing our own less than (lt). ex. I can use it to compare two circles. The default python is for numeric type (5 < 6) or (4.3<4.5)
# __str__(self) is a string representation of this class, if I use print(object), whatever is written inside __str__ method will be printed
# The following class is for split events, a named tuple _SplitEvent is passed into the class
class _SplitEvent(namedtuple("_SplitEvent", "distance, intersection_point, vertex, opposite_edge")):
__slots__ = ()
def __lt__(self, other):
return self.distance < other.distance
def __str__(self):
return "{} Split event @ {} from {} to {}".format(self.distance, self.intersection_point, self.vertex, self.opposite_edge)
# The following class is for Edge events, a named tuple _SplitEvent is passed into the class
class _EdgeEvent(namedtuple("_EdgeEvent", "distance intersection_point vertex_a vertex_b")):
__slots__ = ()
def __lt__(self, other):
return self.distance < other.distance
def __str__(self):
return "{} Edge event @ {} between {} and {}".format(self.distance, self.intersection_point, self.vertex_a, self.vertex_b)
In [15]:
# Two other namedtuples, _OriginalEdge for storing Original edges with left and right bisector. The other one is Subtree to store the skeleton arcs as source (intersection point/node), height (distance from the event edge), and sinks (neighbours/connections of the source node).
_OriginalEdge = namedtuple("_OriginalEdge", "edge bisector_left, bisector_right")
Subtree = namedtuple("Subtree", "source, height, sinks")
If I give the vertices in the counter clockwise (CCW) order as (100, 50), (200, 200), ........,(40,520), (40,50). I get the following polygon
In polyskel github readme.md documentation, the vertices are given in the following order (40, 50), (40, 520), (625,425), (500,325), (635,250), (635,10), (250,40), (200,200), (100,50), which is in clockwise (CW) order. In code comments, it is clearly mentioned the order is counterclockwise. In addition, the original image in polyskel implementation isshown as below:
If you look at them closely, the polyskel documentation image is mirrored in the y-axis. i.e. the positive x axis runs from left of the screen to the right of the screen and the positive y-axis runs from top of the screen to the bottom of the screen. If you assume this convention, then you get the following image which is similar to polyskel documentation. The order of vertices is also in counterclockwise (CCW) direction.
This convention is also confirmed by the fact that the cross product signs are reversed. If you take polyskel convention, the cross product vector of X with Y axis is perpendicular to the plane of paper/screen and the postitive direction for 'z' is into the screen (pointing away from you). In normal convention, the cross product vector of X with Y axis is perpendicular to the plane of paper/screen but the positive direction for 'z' is away from the screen (pointing towards you).
First, let me run the code with polyskel convention to make sure we understand the code. Then we can run it according to normal convention
In [16]:
# One of the most important classes of this program is _LAVertex, I will try to explain this class as detailed as possible
class _LAVertex:
# you pass the point (vertex), left edge for the vertex, right edge. if possible, direction vectors also. You initiate the vertex validity, self._valid = True
def __init__(self, point, edge_left, edge_right, direction_vectors=None):
self.point = point
self.edge_left = edge_left
self.edge_right = edge_right
self.prev = None
self.next = None
self.lav = None
self._valid = True
# TODO this might be handled better. Maybe membership in lav implies validity?
#The order of vertices of the polygon must be given in counter clockwise direction, so an edge pointing towards or entering the vertex is the left edge and an edge pointing away from or leaving the vertex is called the right edge
#Edge left direction is reversed (multiplied by -1) so that the direction is pointing away from or leaving the vertex, by doing this we can see that both edge left and edge right start from the vertex, these vectors are called as creator vectors
#If there are no direction vectors, creator vectors formed from edge left and edge right vectors will be assigned to direction vectors
#creator_vectors = (edge_left.v.normalized() * -1, edge_right.v.normalized())
creator_vectors = (edge_left.v.normalized()*-1, edge_right.v.normalized())
if direction_vectors is None:
direction_vectors = creator_vectors
#To find whether a vertex is reflex or not, we take the cross product of vectors at the vertex (direction vectors). For reflex vertex, the angle between reversed left edge and right edge is less than 180 which makes the cross product resultant vector magnitude positive but the resultant vector direction is perpendicular to screen and points towards you, which according to polyskel convention is negative, so if the cross product of creation/direction vectors is negative, then the vertex is reflex
self._is_reflex = ((_cross(*direction_vectors)) < 0) # original polyskel code convention
#self._is_reflex = ((_cross(*direction_vectors)) > 0) # Vasan/normal code convention
#Angular bisector and description of vector geometry is given in the markdown cell below this cell
self._bisector = Ray2(self.point, operator.add(*creator_vectors) * (-1 if self.is_reflex else 1))
#log.info("Created vertex %s", self.__repr__())
#_debug.line((self.bisector.p.x, self.bisector.p.y, self.bisector.p.x + self.bisector.v.x * 100, self.bisector.p.y + self.bisector.v.y * 100), fill="blue")
#@property decorator helps us to change the attributes (bisector,_is_reflex,.lav._slav._original_edges) and ensures that the change is reflected on other attributes which depend on this attribute. Without property decorator, you can achieve that by making these attributes as part of methods. If you do that, you have to call them as bisector() not self.bisector. If you forget () sometimes, it might break the code. @property decorator is a convenient way. For further explanation, please have a look at this link - https://www.machinelearningplus.com/python/python-property/
@property
def bisector(self):
return self._bisector
@property
def is_reflex(self):
return self._is_reflex
@property
def original_edges(self):
# pylint: disable=maybe-no-member
return self.lav._slav._original_edges
# next_event is an important method which calculates the events for a given a vertex
def next_event(self):
events = []
if self.is_reflex:
# a reflex vertex may generate a split event
# split events happen when a vertex hits an opposite edge, splitting the polygon in two.
#log.debug("looking for split candidates for vertex %s", self)
for edge in self.original_edges: #self.original_edges calls the function original edges inside this _LAVertex class which returns _slav._original_edges
if edge.edge == self.edge_left or edge.edge == self.edge_right:
continue
#log.debug("\tconsidering EDGE %s", edge) #uncomment in polyskel code
# a potential b is at the intersection of between our own bisector and the bisector of the
# angle between the tested edge and any one of our own edges.
# we choose the "less parallel" edge (in order to exclude a potentially parallel edge) #explanation in cells below
leftdot = abs(self.edge_left.v.normalized().dot(edge.edge.v.normalized()))
rightdot = abs(self.edge_right.v.normalized().dot(edge.edge.v.normalized()))
selfedge = self.edge_left if leftdot < rightdot else self.edge_right # if left edge is less parallel, choose it as self edge and right edge as other
otheredge = self.edge_left if leftdot > rightdot else self.edge_right # if right edge is less parallel, choose it as other edge and right edge as self edge
i = Line2(selfedge).intersect(Line2(edge.edge)) #intersection point, i between less parallel edge (self edge) and edge to the checked(edge.edge)
if i is not None and not _approximately_equals(i, self.point): #check that intersection exists(not None) and not equal to the reflex vertex point(self.point)
# locate candidate b
linvec = (self.point - i).normalized() #(vertex point - intersection point) normalized
edvec = edge.edge.v.normalized() # normalized edge vector
if linvec.dot(edvec) < 0: #proper intersection happens only if linvec.dot(edvec) > 0, else, the direciton of edvec is reversed (look at cell below)
edvec = -edvec
bisecvec = edvec + linvec # bisector vector of edvec and linvec
if abs(bisecvec) == 0: # To make sure that linvec and edvec are not parallel
continue
bisector = Line2(i, bisecvec) #Line2 in euclid extends infintely on both sides, bisector passes through i (intersection) and has a direction of bisecvec
b = bisector.intersect(self.bisector) #find the intersection of reflex vector bisector (self.bisector) with bisecvec
if b is None: #if there is no intersection, then move on to the next vertex
continue
# If there is an intersection, we need to check the eligibility of b
# a valid intersection b should lie within the area limited by the edge and the bisectors of its two vertices (look at the note in below cells):
xleft = _cross(edge.bisector_left.v.normalized(), (b - edge.bisector_left.p).normalized()) > 0 #sign is reversed because of polyskel convention
xright = _cross(edge.bisector_right.v.normalized(), (b - edge.bisector_right.p).normalized()) < 0 #sign is reversed - polyskel convention
xedge = _cross(edge.edge.v.normalized(), (b - edge.edge.p).normalized()) < 0 #sign is reversed - polyskel convention
if not (xleft and xright and xedge): #if any of the conditions are not satisfied, then it is invalid, so move on to the next vertex
#log.debug("\t\tDiscarded candidate %s (%s-%s-%s)", b, xleft, xright, xedge)
continue
#log.debug("\t\tFound valid candidate %s", b)
events.append(_SplitEvent(Line2(edge.edge).distance(b), b, self, edge.edge)) #append it to the events using _SplitEvent named tuple
i_prev = self.bisector.intersect(self.prev.bisector) #intersection of current convex vertex bisector with previous vertex bisector
i_next = self.bisector.intersect(self.next.bisector) #intersection of current convex vertex bisector with next vertex bisector
if i_prev is not None: #i_prev exists only for convex vertex
events.append(_EdgeEvent(Line2(self.edge_left).distance(i_prev), i_prev, self.prev, self)) #append the event to events with _EdgeEvent tuple (three items) with distance between i_prev (intersection point) and edge_left line, i_prev (intersection point) and previous vertex
if i_next is not None: #i_next exists only for convex vertex
events.append(_EdgeEvent(Line2(self.edge_right).distance(i_next), i_next, self, self.next)) #append the event to events with _EdgeEvent tuple (three items) with distance between i_next (intersection point) and edge_right line, i_next (intersection point) and next vertex
if not events: #If there are no events, return None
return None
ev = min(events, key=lambda event: self.point.distance(event.intersection_point)) #Find the minimum of the events by distance. The procedure for convex events is given in the cell after bisectors handwritten note
#log.info("Generated new event for %s: %s", self, ev)
return ev
#method to invalidate the vertex, so it is not considered the next time
def invalidate(self):
if self.lav is not None:
#pylint: disable=too-many-function-args
self.lav.invalidate(self)
else:
self._valid = False
@property
def is_valid(self):
return self._valid
def __str__(self):
return "Vertex ({:.2f};{:.2f})".format(self.point.x, self.point.y)
def __repr__(self):
return "Vertex ({}) ({:.2f};{:.2f}), bisector {}, edges {} {}".format("reflex" if self.is_reflex else "convex",
self.point.x, self.point.y, self.bisector,
self.edge_left, self.edge_right)
Convex Event Explanation¶
In [20]:
# The first class in skeletonize function
class _SLAV:
def __init__(self, polygon, holes):
contours = [_normalize_contour(polygon)] #_normalize_contour function explained in previous cells for a polygon
contours.extend([_normalize_contour(hole) for hole in holes]) #Do the _normalize contour function for the number of holes in polygon, add them to the contours list
self._lavs = [_LAV.from_polygon(contour, self) for contour in contours] # Now call the from_polygon() method _LAV class (one cell below this cell)
# store original polygon edges for calculating split events, the _OriginalEdge Tuple stores three items for every vertex - 1.The linesegment from previous vertex to the current vertex 2. Bisector of previous vertex 3. Bisector of current vertex. chain.from_iterable(lists) is a method of itertools module. It combines the various lists to make a single list, so you can iterate through each element of the lists. Here we iterate through every vertex in the self._lavs
self._original_edges = [
_OriginalEdge(LineSegment2(vertex.prev.point, vertex.point), vertex.prev.bisector, vertex.bisector)
for vertex in chain.from_iterable(self._lavs)
]
# An iterator method which iterates through the list of active vertices (LAV)
def __iter__(self):
for lav in self._lavs:
yield lav
# A method for finding the length of _lavs
def __len__(self):
return len(self._lavs)
# Check whether the _lavs is empty or not
def empty(self):
return len(self._lavs) == 0
# Handling of edge events
def handle_edge_event(self, event):
sinks = [] #initialize an empty list called sinks
events = [] #initialize an emptly list called events
lav = event.vertex_a.lav #assign the lav of vertex_a of the event to lav
if event.vertex_a.prev == event.vertex_b.next: #if prev of vertex_a is equal to next of vertex_b then it is a peak event
#log.info("%.2f Peak event at intersection %s from <%s,%s,%s> in %s", event.distance, event.intersection_point, event.vertex_a, event.vertex_b, event.vertex_a.prev, lav)
self._lavs.remove(lav) #remove the lav
for vertex in list(lav): #for every vertex in the list of lav
sinks.append(vertex.point) # append the vertex point to sink
vertex.invalidate() #invalidate the vertex
else:
#log.info("%.2f Edge event at intersection %s from <%s,%s> in %s", event.distance, event.intersection_point, event.vertex_a, event.vertex_b, lav)
new_vertex = lav.unify(event.vertex_a, event.vertex_b, event.intersection_point) #pass vertex_a, vertex_b, intersection_point of event to unify method of LAV class. The unify method returns the replacement vertex which is an instance of _LAVertex with intersection point
if lav.head in (event.vertex_a, event.vertex_b): #check if lav.head is in event.vertex_a or event.vertex_b
lav.head = new_vertex #in that case, the new_vertex (replacement vertex) becomes the head
sinks.extend((event.vertex_a.point, event.vertex_b.point)) #vertex_a.point and vertex_b.point are added to the sinks
next_event = new_vertex.next_event() #now calculate next_event for this new_vertex(replacement vertex)
if next_event is not None: #if there is a next_event for new_vertex (replacement vertex)
events.append(next_event) #append it to the events, which was initially an empty list
return (Subtree(event.intersection_point, event.distance, sinks), events) #return both - a namedtuple Subtree and events
def handle_split_event(self, event):
lav = event.vertex.lav
#log.info("%.2f Split event at intersection %s from vertex %s, for edge %s in %s", event.distance,
#event.intersection_point, event.vertex, event.opposite_edge, lav)
sinks = [event.vertex.point] #for split event, the only sink will be the reflex vertex
vertices = []
x = None # right vertex
y = None # left vertex
norm = event.opposite_edge.v.normalized() #normalized opposite edge of event
for v in chain.from_iterable(self._lavs): #Take all the vertices of polygon
#log.debug("%s in %s", v, v.lav)
#check whether the vertex v's normalized left edge vector is equal to norm and the point of event opposite edge (P) is equal to vertex v's left edge point
#else if whether the vertex v's normalized right edge vector is equal to norm and the point of event opposite edge (P) is equal to vertex v's right edge point
if norm == v.edge_left.v.normalized() and event.opposite_edge.p == v.edge_left.p:
x = v #if true assign the vertex v to x(right vertex)(initially none)
y = x.prev #the previous of x is pointed to y (left vertex) (initially none)
elif norm == v.edge_right.v.normalized() and event.opposite_edge.p == v.edge_right.p:
y = v #elif true assign the vertex v to y(left vertex)
x = y.next #the next of y is pointed to x(right vertex)
if x: #if x is not none, check the direction of cross product vectors
xleft = _cross(y.bisector.v.normalized(), (event.intersection_point - y.point).normalized()) >= 0
xright = _cross(x.bisector.v.normalized(), (event.intersection_point - x.point).normalized()) <= 0
#log.debug("Vertex %s holds edge as %s edge (%s, %s)", v, ("left" if x == v else "right"), xleft, xright)
if xleft and xright: # if both are true exit the iteration loop
break
else: # else assign x and y as None
x = None
y = None
if x is None: #if x is None
#log.info("Failed split event %s (equivalent edge event is expected to follow)", event)
return (None, []) #A failed split event and exits the function
v1 = _LAVertex(event.intersection_point, event.vertex.edge_left, event.opposite_edge) #instance of _LAVertex with event's intersection, edge left and opposite edge. It means that in v1 instance for the intersection point, edge left acts as left edge and opposite edge acts as right edge
v2 = _LAVertex(event.intersection_point, event.opposite_edge, event.vertex.edge_right) #instance of _LAVertex with event's intersection, opposite edge and edge right. It means that in v2 instance for the intersection point, opposite edge acts as left edge and edge right acts as right edge
#Have a look at the handwritten notes to get a better idea
v1.prev = event.vertex.prev # for v1, the previous vertex of event vertex is also the previous vertex of v1
v1.next = x #the next vertex of v1 is the right vertex x
event.vertex.prev.next = v1 # the next vertex of event.vertex.prev is the v1 vertex itelf
x.prev = v1 # v1 vertex is the previous vertex of right vertex x
v2.prev = y # the previous vertex of v2 is the left vertex y
v2.next = event.vertex.next #the next vertex of v2 is the next vertex of event vertex
event.vertex.next.prev = v2 # event.vertex.next's previous vertex is v2
y.next = v2 # the next vertex of y is the vertex v2
new_lavs = None # Initially there are no new_lavs
self._lavs.remove(lav) #remove the current lav
if lav != x.lav: #if event vertex lav is not equal to x.lav (which means both are not same points)
# the split event actually merges two lavs
self._lavs.remove(x.lav) #remove the x.lav
new_lavs = [_LAV.from_chain(v1, self)] #class method from_chain in _LAV class is called without instantiation, which is the advantage of class method
else:
new_lavs = [_LAV.from_chain(v1, self), _LAV.from_chain(v2, self)] #create two new _LAVs from vertex v1 and v2
for l in new_lavs:
#log.debug(l)
if len(l) > 2:
self._lavs.append(l) #append the new lavs to the original ones
vertices.append(l.head) #append the head vertex of each lav
else: # if length of new_lavs is less than 2
#log.info("LAV %s has collapsed into the line %s--%s", l, l.head.point, l.head.next.point)
sinks.append(l.head.next.point) # append the next vertex point of head of lav
for v in list(l): #invalidate the vertex is list l (if you have less than two lavs)
v.invalidate()
events = [] # initialize a null list events
for vertex in vertices: #for vertex in the head
next_event = vertex.next_event() #find next event
if next_event is not None:
events.append(next_event) #append the next event to events list
event.vertex.invalidate() #invalidate the event vertex
return (Subtree(event.intersection_point, event.distance, sinks), events) #return the subtree
